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Many simulations of gravitational collapse to black holes become inaccurate before the total emitted gravitational radiation 
can be determined. The main difficulty is that a significant component of the radiation is still in the near-zone, strong field 
region at the time the simulation breaks down. We show how to calculate the emitted waveform by matching the numerical 
simulation to a perturbation solution when the final state of the system approaches a Schwarzschild black hole. We apply the 
technique to two scenarios: the head-on collision of two black holes, and the collapse of a disk to a black hole. This is the first 
reasonably accurate calculation of the radiation generated from colliding black holes that form from matter collapse. 
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I. INTRODUCTION 

One of the most significant goals of numerical rela- 
tivity is to calculate the gravitational wave forms from 
realistic, astrophysically interesting systems. Until now, 
obtaining a good estimate of the asymptotic wave form 
has required long numerical evolutions so that the emit- 
ted waves can be propagated far from the source. There, 
either radiation extraction techniques can be used, or 
the solution can be matched onto a null-cone integration. 
Unfortunately, long evolutions are difficult for a number 
of reasons. These include throat stretching when black 
holes form, numerical instabilities associated with curvi- 
linear coordinate systems, and the effects of approximate 
outer boundary conditions. Moreover, evolutions using 
certain time-slicing conditions do not allow simulations 
such as collisions to progress far enough for the black 
holes to completely merge; rather they asymptote to a 
finite separation. At this time, no numerical algorithm 
is known that deals effectively with all these problems in 
the context of multidimensional, radiating spacetimes. 

Recently however, an important idea introduced by 
Price and PuUin j^] provides a possible way to circum- 
vent some of these problems in certain circumstances. 
Price and Pullin considered Misner initial data repre- 
senting two black holes at a moment of time-symmetry. 
They realized that when the two holes were sufficiently 
close, the system could be treated as a single, perturbed 
black hole. By applying gauge invariant perturbation 
theory, they calculated this perturbation and evolved it 
using the Zerilli equation. This allowed them to compute 
asymptotic wave forms and emitted energies. Remark- 
ably, for fairly small separations, the energies and wave 
forms agreed well with the results of fully relativistic nu- 
merical simulations Q . In a subsequent study, Abrahams 
and Cook extended this technique and applied it to 
initial data representing boosted black holes with a com- 
mon apparent horizon. The radiated energies computed 
from these initial data sets agreed quite well with those 



from the fully relativistic simulations for cases in which 
the initial separations were large. 

In this paper, we apply the same technique to numer- 
ically evolved spacetimes that at late times can be ap- 
proximated as a single perturbed black hole. The per- 
turbation of the black hole is computed on a spatial slice 
from the numerical metric and extrinsic curvature using 
gauge-invariant perturbation theory. Wave forms are de- 
termined by evolving this perturbation to infinity using 
the Zerilli equation. 

The important difference between this method and 
standard finite-radius radiation extraction techniques (cf. 
[Q) is the use of the horizon as an inner boundary at 
which radiation is purely ingoing. Once the horizon 
forms, and the spacetime settles down to a quasi-static 
state, data from the numerical simulation is no longer re- 
quired. The simulation can be terminated at a relatively 
early epoch; it is not necessary to propagate the radiation 
pulse out to the weak-field regime. Another advantage of 
the perturbation method is that the effects of backscat- 
ter of the waves off the black hole curvature are auto- 
matically taken into account via the Zerilli equation in- 
tegration of the perturbation to large radii. These effects 
can only be incorporated approximately with standard 
extraction methods in which an integration over a time- 
like cylinder is used to separate off near-zone effects. It 
should be emphasized that the perturbation approach re- 
quires the formation of a black hole during the numerical 
simulation. The black hole horizon cuts off the evolution 
inside the black hole from influencing the waves outside. 
The method would be inadequate for, say, an oscillating 
neutron star spacetime, where there is no horizon. 

In Sec. |l| we briefly describe the numerical relativ- 
ity code and detail the implementation of the gauge- 
invariant black hole perturbation method used to com- 



pute wave forms. In Sec. Ill we demonstrate the applica- 



tion of this method to calculations of black hole collisions 
and disk collapse. 
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II. METHODOLOGY 

The numerical spacetimes representing colliding black 
holes and collapsing disks are generated using a code that 
solves the field equations of general relativity coupled to 
a collisionless-particle matter source. Full details of the 
equations solved and the numerical method can be found 
in Ref. |^,^ The code employs the maximal time slic- 
ing condition and the quasi-isotropic spatial gauge, and 
applies to nonrotating spacetimes in axisymmetry. The 
metric, expressed in spherical polar coordinates, takes 
the form 
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The code solves the Hamiltonian constraint to compute 
0g-»7/3^ solves the momentum constraints for the compo- 
nents and fT^ of the extrinsic curvature, and evolves 
the variables 77 and Kg. The kinematical variables a and 

are determined from the maximal slicing condition and 
spatial gauge conditions respectively. In Ref. j|] the code 
was used to study head-on collisions and in Ref. (D it was 
used to study disk collapse. 

A thorough discussion of the perturbed black hole ap- 
proximation will be published elsewhere |Q. Here we 
sketch the general idea and list the equations relevant 
to the current application. The theory of gauge in- 
variant perturbations of Schwarzschild was first devel- 
oped by Moncrief Moncrief showed how to use 
the Regge- Wheeler multipole amplitudes and their ra- 
dial derivatives to construct unconstrained even and odd- 
parity functions that are invariant under infinitesimal co- 
ordinate transformations about the Schwarzschild back- 
ground. These functions satisfy the Zerilli or Regge- 
Wheeler equation for the appropriate ^-mode. For our 
approach, which involves extraction of the perturbation 
data on a spatial slice, the background Schwarzschild 
time-derivative of the gauge invariant function is also re- 
quired in order to pose initial data for integration of the 
second-order Zerilli equation (see (8|j|,01. 

From the viewpoint of our approximation, the non- 
spherical part of the three-metric, the shift, and the 
extrinsic curvature are treated as perturbations about 
a static Schwarzschild black hole in the (approxi- 
mate) exterior region r > M/2, where M is the the 
ADM mass of the slice. In quasi-isotropic gauge, the 
static Schwarzschild metric takes the standard isotropic- 
coordinate form (/3i = 0, 77 = 0, = 1-1- M/2r, 
Kij — 0). To compute the perturbation using standard 
gauge-invariant formalism, first we make the background 
coordinate transformation of the numerically computed 
metric and extrinsic curvature on a given spatial slice 
from isotropic to Schwarzschild radial coordinates r — > i?, 
where R = r(l -I- M/2r)^. Then we compute the the fol- 
lowing Regge- Wheeler multipole amplitudes by integrat- 
ing over coordinate two-spheres at constant R, 
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where x = cos 9. Note that we have specialized to ax- 
isymmetry, so the azimuthal quantum number m — 0. 
We also compute the normal Lie-derivatives of the mul- 
tipole amplitudes: 
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From these we form the gauge-invariant function Qi 
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where iV^ ^ 1 - 2M/r and k = {£ - l){£ + 2) + 6M/R. 
We also compute its Schwarzschild time-derivative Qe = 
NCnQe, where CnQi is constructed by replacing the am- 
plitudes appearing in Eq. ( p^ ) with their normal Lie- 
derivatives. Note that it is valid to compute the time 
derivative using NC^- the difference between these two 
operators is proportional to the Lie-derivative along the 
shift, Cf3^ which is higher order. The two functions Qe 
and Qi provide initial data for integration of the Zerilli 
equation. 



(11) 



where r^, = R + 2M \n{R/2M — 1) is the tortoise coordi- 
nate and the Zerilli potential is given by 



ViiR) = 
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FIG. 1. The radiated energy as a fraction of the total 
mass-energy, determined with the perturbation method, as 
a function of coordinate time during a simulation of colliding 
clusters (Case bl). The curve changes from dashed to solid 
at the time an apparent horizon encompasses both clusters. 
Note that once this common horizon forms, the amount of 
radiation emitted does not depend on when the perturbation 
calculation is performed. 
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Numerically, the initial data [Q i{R) ^ Q ^{R)) is inter- 
polated onto a fine mesh ranging from = —500 to 
r* = 1000 and the Zerilli equation is integrated using a 
2nd-order accurate leapfrog scheme until the perturba- 
tion has all been propagated far away from the peak of 
the potential, \r^\ — > oo. At large radius the function 
approaches the even-parity gravitational wave amplitude 



rhj^ I sin^ Q = 
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III. NUMERICAL RESULTS 

Applying the above equations to numerically generated 
spacetimes is have criteria for when we expect the method 
to give reliable results. Typically, we wait until an appar- 
ent horizon has formed and the distortions of this horizon 
have become fairly small. In the case of collisions, we wait 
until a common apparent horizon, encompassing both ob- 
jects, has formed. We quantify the requirement of small 



distortions by requiring that the proper polar and equato- 
rial circumferences of the apparent horizon be within 10% 
of their Schwarzschild values, i.e., 0.9 ^ Cpo//47rAf < 1.1 
and, 0.9 Ceq/^T^M ^ 1.1. We also require that the area 
of the apparent horizon be within 5% of its Schwarzschild 
value: 0.95 < ^A^/IQttA'P < 1.05. 

Consider a concrete example that illustrates our ap- 
proach. In Fig. II] we show the radiated energy extracted 
when using the perturbation method at different times 
during a numerical simulation. The black hole collision 
shown. Case bl, is discussed below. Fig. Ill A shows the 



particle positions and apparent horizon structure at four 
times during the course of the evolution. The waveform 
extracted using the perturbation method at each of those 
times is also shown. After the common apparent horizon 
forms at a time of about t/M ~ 5.0, the waveform shape 
and amplitude become fairly constant. The disjoint hori- 
zons form earlier, at t/M ~ 4.0. The efficiency - the total 
radiated energy divided by the total mass of the space- 
time - is shown in Fig |l|. This also levels off to a nearly 
constant value once the black holes are encompassed by 
a common horizon. This effect is to be expected if the 
calculation is at all meaningful: the final wave forms and 
energy emission should not depend on when the numeri- 
cal data is extracted for the perturbation equations once 
the final black hole has formed and settled down. Eventu- 
ally, the effects of throat stretching cause the numerical 
solution to become inaccurate and the radiated energy 
begins to diverge rapidly. In all cases where the hori- 
zon settles down to being nearly spherical while the full 
numerical solution is still accurate (judged by the con- 
stancy of the Brill mass |^), the radiated energy remains 
constant as a function of time as required. 

We did encounter some cases where is was impossible 
to determine radiated energy to better than a factor of 
2. Generally this was because the numerical solution be- 
came inaccurate before the horizon settled down. These 
situations are exacerbated by the fact that the quasi- 
isotropic spatial gauge conditions and the Hamiltonian 
and momentum constraints solved in our algorithm are 
elliptic equations. Hence, they instantaneously propa- 
gate numerical errors to the entire solution, even those 
occurring inside the black hole horizon. For this reason, 
even in the relatively accurate cases, we focus on the 
quadrupole, £ = 2 radiation. At present we do not have 
confidence in the lower amplitude i = 4 wave forms as 
the signals are comparable in size to the numerical errors. 

The table lists the 14 cases of colliding clusters and disk 
collapses that were carefully studied using the perturba- 
tion technique. For each case we give the parameters of 
the initial data including proper separations where appli- 
cable. The radiation efficiency is computed for each case. 
In general, we were able to determine the radiated energy 
to better than 20% accuracy, and in many instances to 
better than < 5%. In the remainder of this section we de- 
tail these simulations and discuss the emitted wave forms 
and radiated energies. 
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FIG. 2. The perturbation method apphed at different 
times during a numerical simulation. The lower four frames 
show the particle positions for a collision of initially boosted 
star clusters (Case bl). The dashed lines indicate disjoint 
marginally outer-trapped surfaces, while the solid line shows 
the common apparent horizon. The upper four frames (corre- 
sponding to the same times) show the asymptotic waveform 
computed using the perturbation method as a function of re- 
tarded time in units of M. Radial coordinates are also in 
units of M. 



A. Collisions of clusters and black holes 

An alternative to the vacuum, topological black hole 
approach to simulating black hole collisions is the col- 
lision of black holes that form from collapsing matter. 
Here we consider nonequilbirium spheroids of collision- 
less matter that separately undergo collapse before hit- 
ting each other head-on. An advantage of such an ap- 
proach is that the topology is Euclidean so that simple 
coordinate systems can be employed. It is also possible 
to study different scenarios including those in which the 
clusters of matter do not individually have marginally 
outer-trapped surfaces surrounding them until after they 
are both inside a common event (and possibly apparent) 
horizon. The major disadvantage of this approach is that 
the matter source must be accurately tracked for the en- 
tire simulation - even when it is inside the black hole. 
This can be very computationally intensive. Of course, 
horizon boundary conditions that would allow one to ex- 
cise the black hole interiors from the numerical simulation 
would solve this problem. 

In Sec. H we have briefly described the fully general rel- 
ativistic collisionless matter evolution code used in these 
simulations. Full details of this mean-field, particle simu- 
lation code, as well as initial data and typical simulations 
of colliding clusters can be found in Ref. [|. All simu- 



lations discussed here were performed on a 200 radial 
by 32 angular zone grid with the matter source sampled 
with 3000 particles. As mentioned above, we monitor the 
constancy of quasi-local mass indicators such as the Brill 
mass for a check on accuracy. For the cases shown these 
resources were sufficient to keep the Brill mass constant 
within a few percent. This constancy is crucial for the 
perturbation approach to give consistent results. 

First we consider initial data representing two clus- 
ters of collisionless matter initially at rest. Since this 
initial data is conformally flat and at a moment of time- 
symmetry, it is only necessary to specify the cluster radii 
and separation and solve the Hamiltonian constraint to 
obtain consistent data. The analytic solution for the con- 
formal factor (f) is given in Ref. Interestingly, in the 
vacuum exterior it agrees with the solution given by Mis- 
ner and Wheeler |^ for two black holes in vacuum. This 
black hole solution can be thought of as having a a three- 
sheeted topology with two Einstein-Rosen bridges each 
joining a separate asymptotically-flat sheet to a single 
upper asymptotically- flat sheet. For cases in which each 
cluster initially resides inside its own apparent horizon 
one would expect identical evolutions to result from the 
cluster initial data and the 3-sheeted data. 

Since most numerical studies(cf. Ref. [^) of two black 
hole spacetimes have used two-sheeted Misner data where 
an isometry is imposed between the upper sheet and the 
single lower sheet, it would be interesting to compare the 
radiation expected from evolution of two-sheeted initial 
data with the cluster data to determine the importance 
of topology. Since both types of data represent the same 
physical situation - colliding black holes - it would be 
disturbing if the differences were large. In a future study 
we will examine this issue in more detail both for time- 
symmetric and boosted initial data, which is difficult to 
implement numerically in the 3-sheeted vacuum picture. 

We have studied cluster collisions with a range of sep- 
arations and compactness. In Case cO the cluster centers 
have an initial coordinate separation of only 0.40M and 
there is a common apparent horizon on the initial slice. 
Not surprisingly, this case behaves like a single perturbed 
black hole. In the top frame of Fig. [II A we show the 



waveform from this case computed with the perturba- 
tion method at an early time during the simulation. The 
waveform is very reminiscent of that found by Price and 
PuUin for Misner data in the small separation limit. 
After a fairly short transient, the waveform is dominated 
by the lowest order £ — 2 quasi-normal mode oscillations 
with characteristic wavelength A ~ 16. 8M. 

For comparison, in the same figure we show the wave- 
form from Case c5 in which the clusters initially have 
a coordinate separation of I.IM. Here there are ini- 
tially no common or disjoint apparent horizons. The 
common horizon forms at a coordinate time of about 
2.0M. The horizon is initially very distorted, Ceq — 0.62, 
Cpoi — 1.16. The perturbation calculation of the wave- 
form and radiated energy is not carried out till t ~ 5.3M 
at which time the horizon has become nearly spheri- 
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FIG. 3. Wave forms from Case cO (upper frame) and Case 
c5 (lower frame) are compared. The asymptotic gravitational 
wave amplitude computed with the perturbation method is 
plotted as a function of retarded time in units of M. 



FIG. 4. The gravitational waveform from Case b3 (dashed 
curve) is compared with the waveform from Case b4 (solid 
curve). The asymptotic wave amplitude calculated with the 
perturbation method is plotted as a function of retarded time 
in units of M. 



cal. Not surprisingly, the wave amplitude is much larger 
in this case since the initial separation is considerably 
larger. The waveform shapes are similar but not iden- 
tical because the cluster infall and collapse are entirely 
hidden inside the black hole in Case cO. 

We have also studied initial data, Cases c2 and c4, for 
which there are disjoint apparent horizons surrounding 
each cluster. These cases can be thought of as true black 
hole collisions. Since the clusters are more compact, the 
wave amplitudes in these cases are considerably larger 
than those from cases with the same coordinate separa- 
tion but larger cluster radii. For these cases it is possible 
to compute the initial proper separation of the disjoint 
apparent horizons. It turns out that the proper separa- 
tion is indeed larger for these compact cases. For these 
two cases the radiated energy efficiencies we obtain seem 
to be compatible with the results from head-on collisions 
of time-symmetric topological black holes ^J. Unfortu- 
nately, the throat stretching is very severe in these cases 
and it is difhcult to have confidence in the energies and 
wave forms except for cases with very small initial coor- 
dinate separations. 

We have also considered cases for which the clusters are 
boosted towards each other, like those evolved in Ref. 
Since these data sets are no longer time-symmetric, it is 
necessary to solve both the Hamiltonian and momentum 
constraints for the initial gravitational field. These ini- 
tial value equations are solved by iteration as described 
in Ref. The advantage of the boosted calculations 
is that the clusters can be started at larger coordinate 
separations and still collide before the numerical results 
become unreliable because of the intrinsic collapse of each 



cluster. As can be seen in the table, the boosted Cases 
bl-b4 produce considerably more radiation than the un- 
boosted clusters. With boosts, the empirical limit on 
radiation efficiency of ^ 10~^ for black hole collisions 
resulting from time-symmetric data (cf. Ref. [^) can 
be exceeded. The largest amount of radiation, efficiency 
5 X 10"'^, is produced in Case b2, which has the largest 
separa tion a nd boost velocity. 

Fig. Ill A compares wave forms from two boosted clus- 
ters. Cases b3 and b4. The coordinate separation and 
initial radii of the clusters is the same in each case, only 
the boost velocity is different. Clearly, the wave ampli- 
tude is much greater in the case of the more violent col- 
lision. Case b4. The wave forms from the boosted cases 
all have the same basic qualitative appearance. The ini- 
tial hump, not present in analysis of boosted Misner-like 
data 1^ , probably arises because the clusters are not ini- 
tially black holes. Quasi-normal modes dominate after 
about 20M in all cases. 



B. Disk collapse 

Collapsing disks of coUisionless matter provide an ex- 
cellent test problem for axisymmetric numerical relativity 
codes. Since the source is all in the equatorial plane, the 
matter evolution equations are all one-dimensional. The 
source is felt by the two-dimensional gravitational field 
via jump conditions. When the disk matter is initially 
at rest, the situation provides an interesting analogy to 
Oppcnheimer-Snyder collapse to a black hole in spherical 
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FIG. 5. Particle positions in the equatorial plane are plot- 
ted at several times during a disk collapse, Case d3. The 
apparent horizon is indicated by the solid line in the third 
and fourth frames. 



symmetry, but with the added important feature of grav- 
itational radiation production. Since the gravitational 
field is dynamical in disk collapse, the full machinery of 
numerical relativity is required to follow the evolution, 
while spherical Oppenheimer-Snyder collapse is analytic. 
In Ref. Q the basic equations for this model were given 
and a number of test simulations discussed. Here we 
restrict our attention to cold matter cases that collapse 
to black holes. By cold matter we mean cases in which 
the particles in the disk have no velocity dispersion at 
t — 0. All simulations shown here were performed on 
a 300 radial by 16 angular zone grid with the matter 
source represented by 12000 particles. Even though disk 
collapse is in a sense the simplest radiating system, in no 
case are we able to track the evolution long enough to 
read off the full waveform directly from the simulation. 
Only with the perturbation method can we extract the 
entire waveform. 

We consider three disk collapses. Cases dl-d3 in the 
table. As might be expected the larger the initial radius, 
the more radiation is produced, although the efficie ncies 
are very small, less than 0.1% in all cases. Fig. [II B 
shows a typical collapse. Case d3. The matter starts at 
rest with initial disk radius R/M = 2.0. The locations 



of the particles in the equatorial plane are shown. After 
the matter passes inside a radius of about r/AI ~ 0.5 an 
apparent horizon forms. As in Newtonian gravitation, 
cold disks are unstable to ring formation. As Ref. ||] 
shows, ring formation here mostly occurs after the black 
hole forms. The waveform computed using the pertur- 
bation appro ach at a time of about t/M = 8 is shown 
in Fig. Ill B , A least squares fit to the two lowest-order 



FIG. 6. The gravitational waveform is plotted as a function 
of retarded time in units of M for disk collapse Case d3 (solid 
line). The least-squares fit to the two lowest-order £ — 2 
Schwarzschild quasi-normal modes is also shown (dashed line) 
for t-r,> 17M. 



£ = 2 Schwarzschild quasi-normal modes is also shown. 
After an initial transient feature, lasting for about 15M, 
the fit is excellent in both amplitude and frequency. 
In Fig. Ill B we compare the waveform computed with 



the perturbation technique (using simulation data at a 
time of about t = 4.7M) with that extracted using the 
usual technique Q for Case d2 (see Ref. Q). Note that 
the simulation terminated within about 15A/ after black 
hole formation, leaving little time for waves to reach the 
numerical detector at r = 8M. Only the radiation from 
the infall phase and the first oscillation of the black hole 
are obtainable using the standard method before the sim- 
ulation goes bad. The agreement between the two meth- 
ods up to this time is reasonable considering that no spe- 
cial care has been taken in separating out the near-zone 
field for the infall part of the waveform. This is a delicate 
matter requiring special determination of initial values 
for the spacelike extraction techniques [0, but is au- 
tomatically taken care of in the perturbation approach 
provided that the waveform is read off the Zerilli inte- 
gration at a sufficiently large radius. 



IV. DISCUSSION 

In summary, we have found that the perturbation ap- 
proach to calculating asymptotic wave forms and radi- 
ated energies is effective not only in analyzing initial data 
sets, but can also be used in conjunction with numerical 
relativity simulations. For many scenarios it is difficult to 
evolve black holes long enough to compute accurately the 
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FIG. 7. The gravitational waveform computed using the 
perturbation method for disk collapse Case d2 is plotted as a 

function of retarded time in units of A/. For comparison, the 
gravitational waveform computed at a radius r/M = 8 using 
the standard spacelike radiation extraction technique is also 
shown (points). 
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wave forms produced by their formation and ring-down. 
As long as the simulation can be prolonged to the point 
that a fairly spherical black hole horizon can be identi- 
fied, the perturbation method presented here should give 
reliable results. The method holds great promise for more 
general systems, since it should be possible to extend the 
method to perturbations of a rotating black hole. 
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TABLE I. Radiated energy efficiencies from cluster collisions and cold disk collapse. 



Case 


description 


zo/M 


Ro/M " 


L/M " 


Lah/M " 


v/c " 




cO 




collision 


0.20 


0.15 


6.97 




0.0 


1.5 X 10"^ 


cl 




collision 


0.40 


0.3 


5.1 




0.0 


1.5 X 10"^ 


c2 




collision 


0.40 


0.1 


9.6 


3.1 


0.0 


1.4 X 10"* 


c3 




collision 


0.475 


0.3 


5.7 




0.0 


5.5 X 10"^ 


c4 




collision 


0.50 


0.1 


10.2 


3.6 


0.0 


3. X 10"" 


c5 




collision 


0.55 


0.3 


6.1 




0.0 


7.5 X 10"^ 


c6 




collision 


0.65 


0.3 


6.4 




0.0 


9.0 X 10"^ 


bl 




boosted collision 


1.0 


0.5 


6.2 




0.05 


1.2 X 10"^ 


b2 




boosted collision 


1.4 


0.5 


6.5 




0.15 


5. X 10"^ 


b3 




boosted collision 


0.75 


0.5 


5.3 




0.05 


7. X 10"* 


b4 




boostcxi collision 


0.75 


0.5 


5.1 




0.10 


4. X 10"^ 


dl 




disk 




1.00 








1.5 X 10"* 


d2 




disk 




1.50 








3.3 X 10"* 


d3 




disk 
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''Coordinate displacement up the axis of the cluster center. 

''Coordinate radius of cluster or disk. 

'^Proper separation of cluster centers. 

''Proper separation of disjoint apparent horizons. 

''Boost velocity of cluster as measured by a normal observer. 

'The radiation efficiency. 
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